Robustness Leads Close to the Edge of Chaos in Gene Regulatory Networks 
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Transcriptional dynamics of gene regulatory networks are regulated in highly precise manner, 
despite a fluctuating environment and mutations. We model these dynamics as those of a coupled 
logistic map on a network and design systems which are robust against phenotypic perturbations 
(perturbations in dynamics), as well as systems which are robust against mutation (perturbations 
in network structure). To achieve such a design, we apply a multicanonical Monte Carlo. Analysis 
based on the maximum Lyapunov exponent and parameter sensitivity shows that systems with 
marginal stability, which are regarded as systems at the edge of chaos, emerge when robustness 
against genotypic perturbations is required. This emergence of the edge of chaos is a self-organization 
phenomenon and does not need a fine tuning of parameters. 
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Gene expression in living cells is a highly complex pro- 
cess, achieved through a network of mutual regulation 
between numerous genes. Transcription-translation reg- 
ulation is robust against both phenotypic and genotypic 
perturbations [H-Q • Put differently, gene expression pat- 
terns are kept stable in spite of intrinsic and extrin- 
sic noise (i.e., phenotypic perturbations) and mutations 
(i.e., genotypic perturbations). It seems reasonable to 
think that robustness against phenotypic perturbations 
has been evolutionarily developed for adapting to noisy 
environments. There are also several advantages in hav- 
ing mutational robustness - it buffers against deleterious 
mutations, and it has been recently suggested that mu- 
tational robustness enhances evolvability 0, [H, H| . 

Robustness, as described above, is a fundamental prop- 
erty of life. This point of view naturally gives rise to the 
question: what kind of system emerges when only robust- 
ness is required ? In this paper, we show that systems at 
the edge of chaos are selected with only the requirement 
of mutational robustness. 

It has long been hypothesized that living systems fa- 
vor the edge of chaos, where stability and chaoticity co- 
exist. Originally, Kauffman Q introduced the Boolean 
network model (N-K model) as a model of a gene regu- 
latory network, and proposed the hypothesis that living 
systems prefer the edge of chaos because it allows sys- 
tems to have complex behaviors || . Here we propose an 
alternative scenario, specifically that the requirement of 
having mutational robustness drives living systems to the 
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edge of chaos, regardless of whether or not living systems 
prefer it. 

We propose a coupled map system as an abstract 
model of a gene regulatory network. Unlike the N-K 
model, each element in this model has its own dynamics. 
Assuming that x\ € (— 1, 1) is the gene expression of the 
z-th gene at time step t, the single gene dynamics are 
written as = G(x\). These dynamics mimic mul- 

tiple processes in an expression of a single gene. In the 
presence of N genes, the dynamics of x\ are expressed as 



N 



' + 1 



= (l-e)G(xD + eJ2W ij G(x t j ), 



(1) 



where e is a coupling constant. Wij describes the in- 
teraction strength from gene j to gene i; Wij <G [0, 1] 
satisfies both conditions Wu = and Ylf Wij = 1 
for each i. We call the matrix W, whose ij element 
Here we choose the logistic map 



is Wij, a network 



g(x,a) = 1 



as G(-). We use the model parameters 
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a and e as (a, e) = (1.8, 0.1). This choice indicates that a 
single disconnected gene exhibits chaotic dynamics. We 
impose an additional constraint on W: the number of 
input links k to each gene is fixed. We note that in the 
case of Wij = 1/(N — 1) for all the system becomes 

the globally coupled map (GCM) [y] and it shows highly 
chaotic behaviors at the parameters (a, e) = (1.8, 0.1). In 
contrast to the N-K model, variables of this model take 
continuous values and a linear stability analysis can be 
applied. 

In this model, a network W is regarded as a geno- 
type while an attractor of dynamics x is regarded as 
a phenotype. Our goal is to design networks under the 
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two different design principles: robustness against pheno- 
typic perturbations (i.e., perturbations in the dynamics 
of gene cc) and robustness against genotypic perturba- 
tions (i.e., perturbations on network W). In both cases, 
only network W is tuned. 

Let us start with the first design principle, namely ro- 
bustness against perturbations in dynamics. In other 
words, we are aiming to design a system with a stable 
attractor. For this end, we use Lyapunov exponent anal- 
ysis and a multicanonical Monte Carlo @, Hf ■ 

Once a network W is given, the finite time maximum 
Lyapunov exponent Ai [9( is calculated for the dynami- 
cal system in Eq. ([T]), starting from a given initial state 
x°*pi|. We perform the simulation up to T = 1500 and 
regard the first T" = 1000 steps as transient and discard 
them. 

Multicanonical Monte Carlo @, [H| allows us to sam- 
ple networks with negative Ai, which indicates that a 
network is stable, and to estimate the probability of a 
stable network. Our multicanonical Monte Carlo strat- 
egy adopted in this study is to perform random walks in 
Ai space and thereby to sample rare networks that have 
negative Ai. We define the density of Ai as 



D(X) 



S(Xi(W) 



X) P {W) ILijdWij, 



where 5 is the Dirac 5-function, and p(W) is the proba- 
bility density that a network W appears under random 
sampling. We consider here a network ensemble in which 
p{W) is given by 

p(W) oc IL<5(£ 3 W l3 - 6(W l3 ) - k)6(Wu), (2) 

where 6(z) is a function that satisfies 9{z) = 1 for z ^ 
and 9{z) = for z = 0. A key quantity of multicanonical 
Monte Carlo is the weight function w(X\) of Ai. In this 
study, the Wang and Landau algorithm [l2j is used to 
tune the weight function iu(Ai). Performing Monte Carlo 
sampling with the weight w(Xi), we modify w(\±) step 
by step until the histogram h(X±) of Ai is sufficiently 
flat. In practice, we calculate the density D(Xi) only in 
the prescribed interval A a < Ai < A&. Details of the 
implementation are given in the supplementary material 
(ST1). 

Figure [T] shows the calculated densities -D(Ai) of Ai 
for the fixed input degree k = 2. We also show -D(Ai) 
for k = 3 <~ 5 in the supplementary material (SF1 - 
SF3). Using density D(Xi) in Fig.[TJ the probability that 
networks with negative Ai are observed under random 
sampling is calculated by F(Ai < 0) = J. D(Xi)dX\. We 
estimate P(Ai < 0) with k = 2 ~ 5 and k = N — l, which 
are shown in FigJ5] Each P(Ai < 0) shows that a stable 
attractor becomes inceasingly rare as N or k increases, 
indicating that these systems are in the chaotic phase 
(we define that a system is in the chaotic phase when 
only positive values of Ai appear as N — > oo). These 
results are consistent with the behavior of GCM with 
(a,e) = (1.8,0.1) @. 
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FIG. 1: (COLOR ONLINE)Density of finite time 
Lyapunov exponent under random sampling of networks 
with input degree k = 2. Calculated in the prescribed 
interval -0.28 < Ai < 0.39 for N = 4 ~ 10 and 
-0.1 < Ai < 0.39 for N = 20 and 30. 




network size N 

FIG. 2: (COLOR ONLINE) Network size dependence of 
the probability P(Ai < 0) of networks with negative Ai. The 
logarithms of P(Ai < 0) for k — 2 — 5 decrease linearly or 
slightly faster than linear as functions of N. The logarithm 
of -P(Ai) for k = N — 1 decreases quadratically. 



Using the second design principle, we design systems 
that are robust against genotypic perturbations (i.e., net- 
work perturbations). In other words, we design, using 
multicanonical Monte Carlo, networks W whose trajec- 
tory on the attractor hardly changes when a small net- 
work perturbation 6W is added. We define the param- 
eter sensitivity and use it as a guiding function of the 
robustness. 

Sensitivity analysis using parameter sensitivity has 
been developed and applied in various fields [lj, fl3 |. 
While most of these studies have dealt with continuous 
time systems, we define the parameter sensitivity for dis- 
crete time systems as follows. 

Let us denote the set of elements in W in Eq. ([1} by 
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FIG. 3: (COLOR ONLINE) The two-dimensional density 

D(Ai, 7) of Ai and 7 for N = 10, k = 2 (left) and k = 5 
(right). A large fraction of networks with small 7 degenerate 
in the region indicated by the red rectangles. The black lines 
indicate A = 0. It should be noted that 7 of Ai > diverges, 
when (T m ax — T') — >• 00. 



a vector W. When a small network perturbation SW is 
introduced into W at t = T', the displacement between 
unperturbed trajectory x t (W) and perturbed trajectory 
x t (W + SW) is approximated by (dx' ' / 8W)SW , where 
dx t /dW is a N x Af 2 matrix. We call this matrix the 
and the time evolution of A* is 



sensitivity matrix A*, 
given by 

A t+i 



dF 

dx 



A' 



dF 

dW' 



where dF/dx is the Jacobian matrix and BF/dW is 
the parametric Jacobian matrix. It should be noticed 
that the two trajectories x t {W) and x t (W + SW) co- 
incide for t < T', and thus A* = for t < T'. The 
growth rate of the displacement between x t (W + SW) 
and x t {W) with respect to the perturbation vector SW 
is obtained by A t SW, where SW = SW/\5W\. The 
maximum value of A f <5T^ at time step t is given by the 
maximum singular value a\ of A* matrix. a\ can be ob- 
tained by performing the singular value decomposition of 
A*. <7* diverges for t — > 00 when the maximum Lyapunov 
exponent of the trajectory is positive. On the other hand, 
a\ oscillates or converges to a constant value when the 
maximum Lyapunov exponent is negative. Note that no 
parameters except for W are perturbed in this study. 
Once a network W is given, A* and its a\ are estimated 
for each time step. We define parameter sensitivity 7 as 
the logarithm of an average of a\ along the trajectory: 



T max -1 



7 = In 



Here, we regard the first T" — 1 steps of the trajectory 
as transient, and discard them. We use T 1 — 1000 and 



Our goal is to sample networks with small 7. In order 
to obtain such networks, we perform random walks in Ai 
space with the fixed u>(Ai) estimated above. These ran- 
dom walks facilitate efficient sampling of the networks 
with small 7. We also obtain the two-dimensional den- 
sity .D(Ai,7) °f ^1 ancl 7 as follows: we construct the 
two-dimensional histogram h(\\,^f) through the random 
walks, and, after /i(Ai,7) i s constructed, D(X\,'y) is cal- 
culated by 

D{Xi,l) °c — 77—-. 

iu(Ai) 

Figure |3] shows £)(Ai,7) for N — 10 with k = 2 and 
5. These results show that although networks with small 
7 can take various values for Ax, the vast majority of 
such networks have negative but near zero Ai (see red 
rectangles in Fig. [3]). 

This indicates that, when robustness against network 
perturbations is optimized, networks with negative but 
near zero Ai will appear with high probability. This ap- 
pearance of systems with marginal stability can be in- 
terpreted as self-organization of the edge of chaos. In 
this scenario, the system automatically comes close to 
the edge of chaos without tuning parameters. 

In order to confirm that optimization of mutational ro- 
bustness leads to the emergence of the edge of chaos, we 
perform simulated annealing. Here, we define a parame- 
ter sensitivity T without the linear approximation: 



t=T> 



\x*(W + 8W)-x t (W)\\ 

I sW 



T„ 



1500. 



where () s w represents an average over realizations of 
SW . We minimize this T using the simulated annealing. 
The average is taken over 100 samples of SW, and the 
parameters T' = 1000 and T max — 1500 are used. In each 
step of the simulated annealing, a transition from the 
current state W to a proposed candidate W' is accepted 
if and only if the ratio exp[—/3(T(W')~T(W))] is smaller 
than a random number uniformly distributed in (0,1]. 
Here, the inverse temperature /3 is a increasing function 
of simulation step n. We choose this /3 as j3 = 10n/3 (for 
n < 30000) and (3 = 00 (for n > 30000). Note that the 
function T that we aim to minimize fluctuates due to the 
finite sample size of SW, and thus occasionally a worse 
network could be accepted or a suitable network could 
be rejected, even for f) = 00. 

In Fig|4l we plot Ai and V for networks that are sam- 
pled during the simulated annealing. These results indi- 
cate that most networks obtained in the last half of the 
simulations (n > 30000) are in the region —0.05 < Ai < 
0, and we regard this region as the edge of chaos. 

In summary, using multicanonical Monte Carlo, we 
have shown that systems at the edge of chaos emerge as a 
self-organization phenomenon with only the requirement 
of mutational robustness. We have also performed simu- 
lated annealing and confirmed this scenario. We empha- 
size that no fine tuning of other parameters, such as num- 
ber of input links k or model parameters (a, e), is needed. 
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FIG. 4: (COLOR ONLINE) Time series of sensitivity T 
(upper panel) and the maximum Lyapunov exponent Ai 
(lower panel) in the course of simulated annealing. The red 

and blue lines indicate results for N = Wwithk — 2 and 
TV = lOwithk = 5, respectively. Ai and F are calculated in 
the same simulations. The inverse temperature is set to be 
/3 = oo in the last half of these simulations. 



The emergence of the edge of chaos with the requirement 
of mutational robustness is somehow counterintuitive be- 
cause mutational robustness seems to be positively cor- 
related with dynamical stability. The mechanism of the 
emergence of the edge of chaos, revealed by multicanon- 
ical Monte Carlo, is as follows: when mutational robust- 
ness is required, selected systems need to have Ai < 
because 7 for Ai > diverges as t — > 00. The density 
-D(Ai) is an increasing function for Ai < 0. Therefore, 
the density of networks becomes largest at Ai ~ under 
the condition Ai < (see FigQ]): the red rectangles in 
Figj3] indicate the degeneracy of a large numbers of net- 



works. Due to this degeneracy, systems at the edge of 
chaos are realized in high probability. 

Similar results have been found in recent studies 
; systems that have the ability to reach a stable fixed 
point with transient chaotic behavior appear with only 
the requirement of robustness against genotypic pertur- 
bations. These results can be also interpreted as the 
emergence of the edge of chaos. However, it has not until 
now been explained why such systems are selected. In 
this paper, we have proposed a mechanism for the emer- 
gence of the edge of chaos, namely that the vast major- 
ity of stable networks have marginal stability, and thus 
networks at thcedge of chaos are selected under random 
sampling. It is reasonable to think that this degeneracy 
of marginally stable networks appears whenever parame- 
ters are set in the chaotic phase, in which chaotic systems 
are obtained under random construction of systems for 
TV — > 00. Based on the fact that similar results were 
found in the previous studies most of what we 

discussed here seems not to depend on the details of a 
model. 

Mutational robustness seems to be a natural require- 
ment for living systems. Therefore our scenario is ex- 
pected to apply to real living systems. In fact, several 
recent studies have suggested that gene networks of real 
organisms stay at theedge of chaos (l8l - [2^ |. Our study 
provides an explanation: the requirement of mutational 
robustness drives such organisms to the edge of chaos, 
whether or not staying in such a regime is preferable for 
living systems. 
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